      SUBROUTINE CLOUDW &
     & (qc,u,v,wsig, &
     &  wrk,difunqc,cond,auto,col,pp,pp10,pp01,qcbx,qcby,qccc, &
     &  hk1, &
     &  dtxs,dtys,dxys,ds0,ds12,taudra,hdamp, &
     &  dxdt,dydt,dtxy,dt,dtl,dx2,dy2,filt,vdampt, &
     
     &  nx,nx1,ny,ny1,ns,ns1,nx2,ny2, &
     &  nx1qc,ny1qc,ns1qc,nx1qv,ny1qv,ns1qv, &
     &  nstep,idrmax,numsmo,nchn1,iobptx,iobpty, &
     
     &  tfct,raylei,dohsmo,dovsmt,ptbout,prt)
     
!c-----------------------------------------------------------------------
!c
!c     time integration of the cloud water equation
!c
!c uses wk1
!c-----------------------------------------------------------------------
!c
      use MPI_VARIABLES, only : &
      & nxs2i, nxs2e, nys2i, nys2e

      implicit none

!     Input/output variables

!     Reals
      
      real(8), intent(inout) :: qc(nxs2i:nxs2e,nys2i:nys2e,0:ns1,1:3)
      real(8), intent(inout) :: u(nxs2i:nxs2e,nys2i:nys2e,0:ns1,2:3)
      real(8), intent(inout) :: v(nxs2i:nxs2e,nys2i:nys2e,0:ns1,2:3)
      real(8), intent(inout) :: wsig(nxs2i:nxs2e,nys2i:nys2e,0:ns1,2:3)
      
      real(8), intent(inout) :: wrk(nxs2i:nxs2e,nys2i:nys2e,0:ns1)
      real(8), intent(inout) :: difunqc(nxs2i:nxs2e,nys2i:nys2e,0:ns1)
      real(8), intent(inout) :: cond(nxs2i:nxs2e,nys2i:nys2e,0:ns1)
      real(8), intent(inout) :: auto(nxs2i:nxs2e,nys2i:nys2e,0:ns1)
      real(8), intent(inout) :: col(nxs2i:nxs2e,nys2i:nys2e,0:ns1)
      real(8), intent(inout) :: pp(nxs2i:nxs2e,nys2i:nys2e,1:4)
      real(8), intent(inout) :: pp10(nxs2i:nxs2e,nys2i:nys2e,1:4)
      real(8), intent(inout) :: pp01(nxs2i:nxs2e,nys2i:nys2e,1:4)
      real(8), intent(inout) :: qcbx(1:2,nys2i:nys2e,0:ns1)
      real(8), intent(inout) :: qcby(nxs2i:nxs2e,1:2,0:ns1)
      real(8), intent(inout) :: qccc(1:2,1:2,0:ns1)

      real(8), intent(inout) :: hk1(nxs2i:nxs2e,nys2i:nys2e)
      
      real(8), intent(in) :: dtxs(0:ns1)
      real(8), intent(in) :: dtys(0:ns1)
      real(8), intent(in) :: dxys(0:ns1)
      real(8), intent(in) :: ds0(0:ns1)
      real(8), intent(in) :: ds12(0:ns1)
      real(8), intent(in) :: taudra(0:ns)
      real(8), intent(in) :: hdamp(0:ns1)
      
      real(8), intent(in) :: dxdt
      real(8), intent(in) :: dydt
      real(8), intent(in) :: dtxy
      real(8), intent(in) :: dt
      real(8), intent(in) :: dtl
      real(8), intent(in) :: dx2
      real(8), intent(in) :: dy2
      real(8), intent(in) :: filt
      real(8), intent(in) :: vdampt

!     Integers
      integer(4), intent(in) :: nx
      integer(4), intent(in) :: nx1
      integer(4), intent(in) :: ny
      integer(4), intent(in) :: ny1
      integer(4), intent(in) :: ns
      integer(4), intent(in) :: ns1
      integer(4), intent(in) :: nx2
      integer(4), intent(in) :: ny2
      integer(4), intent(in) :: nx1qc
      integer(4), intent(in) :: ny1qc
      integer(4), intent(in) :: ns1qc
      integer(4), intent(in) :: nx1qv
      integer(4), intent(in) :: ny1qv
      integer(4), intent(in) :: ns1qv      
      integer(4), intent(in) :: nstep
      integer(4), intent(in) :: idrmax
      integer(4), intent(in) :: numsmo
      integer(4), intent(in) :: nchn1
      integer(4), intent(in) :: iobptx
      integer(4), intent(in) :: iobpty

!     Characters

!     Logicals
      logical, intent(in) :: tfct
      logical, intent(in) :: raylei
      logical, intent(in) :: dohsmo
      logical, intent(in) :: dovsmt
      logical, intent(in) :: ptbout
      logical, intent(in) :: prt      

!     Local variables
!     Reals
      real(8), allocatable :: work1(:,:,:)
      real(8), allocatable :: work2(:,:,:)
      real(8), allocatable :: work3(:,:,:)
      real(8), allocatable :: work4(:,:,:)
      real(8), allocatable :: work5(:,:,:)
      real(8), allocatable :: work6(:,:,:)
      real(8), allocatable :: work7(:,:,:)
      real(8), allocatable :: work8(:,:,:)

      real(8) :: tavx
      real(8) :: tavy
      real(8) :: tavs

!     Integers
      integer(4) :: ix ! Loop index
      integer(4) :: iy ! Loop index
      integer(4) :: is ! Loop index
      integer(4) :: k ! Loop index      

      if(tfct) then
         do k=1,3
         do is=0,ns
         do iy=1,ny1
         do ix=1,nx1
            qc(ix,iy,is,k)=qc(ix,iy,is,k)*pp(ix,iy,k)
         enddo
         enddo
         enddo
         enddo
         allocate(work1(0:nx1,0:ny1,0:ns1))
         allocate(work2(0:nx1,0:ny1,0:ns1))
         allocate(work3(0:nx1,0:ny1,0:ns1))
         allocate(work4(0:nx1,0:ny1,0:ns1))
         allocate(work5(0:nx1,0:ny1,0:ns1))
         allocate(work6(0:nx1,0:ny1,0:ns1))
         allocate(work7(0:nx1,0:ny1,0:ns1))
         allocate(work8(0:nx1,0:ny1,0:ns1))
!         call fct3d(qc(0,0,0,3),qc(0,0,0,2),qc(0,0,0,1) &
!     &      ,u(0,0,0,3),u(0,0,0,2) &
!     &      ,v(0,0,0,3),v(0,0,0,2) &
!     &      ,wsig(0,0,0,3),wsig(0,0,0,2) &
!     &      ,nx,ny,ns &
!     &      ,dtxs,dtys,dtxy,dxys,dxdt,dydt,dt,ds0 &
!     &      ,work1,work2,work3,work4,work5,work6,work7,work8)
         deallocate(work1)
         deallocate(work2)
         deallocate(work3)
         deallocate(work4)
         deallocate(work5)
         deallocate(work6)
         deallocate(work7)
         deallocate(work8)

         do k=1,3
         do is=0,ns
         do iy=1,ny1
         do ix=1,nx1
           qc(ix,iy,is,k)=qc(ix,iy,is,k)/pp(ix,iy,k)
         enddo
         enddo
         enddo
         enddo
         do is=1,ns-1
         do iy=2,ny
         do ix=2,nx
           qc(ix,iy,is,3)=qc(ix,iy,is,3)+dtl*(difunqc(ix,iy,is) &
     &       +cond(ix,iy,is)-auto(ix,iy,is)-col(ix,iy,is))/pp(ix,iy,3)
         enddo
         enddo
         enddo
      else
        do is=1,ns-1
        do iy=2,ny
        do ix=2,nx
        tavx=(u(ix,iy,is,2)*pp10(ix,iy,2)*(qc(ix+1,iy,is,2) &
        &  +qc(ix,iy,is,2)) &
        &  -u(ix-1,iy,is,2)*pp10(ix-1,iy,2)*(qc(ix,iy,is,2) &
        &  +qc(ix-1,iy,is,2)))/dx2 
        tavy=(v(ix,iy,is,2)*pp01(ix,iy,2)*(qc(ix,iy+1,is,2) &
        &  +qc(ix,iy,is,2)) &
        &  -v(ix,iy-1,is,2)*pp01(ix,iy-1,2)*(qc(ix,iy,is,2) &
        &  +qc(ix,iy-1,is,2)))/dy2
        tavs=(wsig(ix,iy,is+1,2)*(qc(ix,iy,is+1,2)+qc(ix,iy,is,2)) &
        &  -wsig(ix,iy,is,2)*(qc(ix,iy,is,2)+qc(ix,iy,is-1,2))) &
        &  /ds12(is)
        qc(ix,iy,is,3)=(qc(ix,iy,is,1)*pp(ix,iy,1) &
        &  -dtl*(tavx+tavy+pp(ix,iy,2) &
        &  *(tavs-difunqc(ix,iy,is)-cond(ix,iy,is)+auto(ix,iy,is) &
        &  +col(ix,iy,is))))/pp(ix,iy,3)
        enddo
        enddo
        enddo
      endif


!c surface cloud water flux (bulk parametrization):
!c ?????
!c
!c horizontal boundary conditions:

      call ptbc(nx,ny,ns,qc(0,0,0,3),qcbx,qcby,qccc,iobptx,iobpty &
     &  ,ptbout,prt,nchn1)

      if(raylei) then
        do is=1,idrmax
        do iy=1,ny1
        do ix=1,nx1
          qc(ix,iy,is,3)=qc(ix,iy,is,3)-dtl/taudra(is)*qc(ix,iy,is,2)
        enddo
        enddo
        enddo
      endif

      do iy=1,ny1
      do ix=1,nx1
        qc(ix,iy,ns,3)=qc(ix,iy,ns-1,3)
        qc(ix,iy,0,3)=qc(ix,iy,1,3)
      enddo
      enddo

      if(dohsmo .and. mod(nstep,numsmo).eq.0) then
         call hsmoot(qc(0,0,0,3),hk1,nx1,ny1,ns,1,nx1,1,ny1,1,ns-1 &
     &      ,hdamp)
      endif
      if(dovsmt .and. mod(nstep,numsmo).eq.0) then
         call vsmoot(qc(0,0,0,3),wrk,nx1,ny1,ns,1,nx1,1,ny1,1,ns-1 &
     &      ,vdampt)
      endif

      if(.not.tfct) then
         call aselin(qc(0,0,0,3),qc(0,0,0,2),qc(0,0,0,1) &
     &   ,0,nx1,0,ny1,0,ns,filt)
      endif

      call radbch(qccc,qcbx,qcby,qc(0,0,0,3),qc(0,0,0,2),qc(0,0,0,1) &
     &   ,nx1,ny1,ns1,1,ns-1,1,nx1,1,ny1)

      return
      END SUBROUTINE CLOUDW
